      subroutine nmrvf8(ln,ecmn,ucentr,cntfug,gridx,nx,
     >   jdouble,njdouble,s1,s2,reg,jstart,jstop)
c
      implicit real*8 (a-h, o-z)
c
      dimension jdouble(njdouble),ucentr(nx),reg(nx),gridx(nx),
     >   cntfug(nx)

      x = gridx(jstart)
      dx =  gridx(jstart+1) - x
      f2 = ucentr(jstart)+cntfug(jstart)-ecmn
      h1 = dx * dx
      h2 = h1/12d0

      if (jstart.eq.1) then
C  We get here if the integration is going to start from the first X point.
C  This means that S1 is the solution of the differential equation at X=0
C  S2 is the solution at the first X point. 
         s1 = 0d0
         t1 = 0d0
C  LN = 1 is a special case
         if (ln.eq.1) t1 = -h1/18d0
         if (ecmn.ne.0d0) t1 = t1 * ecmn
      else
         j = jstart-1
         f1 = ucentr(j)+cntfug(j)-ecmn
         t1 = s1 * (1d0-h2 * f1)
      end if

      reg(jstart) = s2
      t2 = s2 * (1d0-h2 * f2)
      
      istart = 2
      do while (jstart.gt.jdouble(istart).and.istart.lt.njdouble)
         istart = istart+1
      end do
      istart = istart-1
      istop = njdouble-1
C  JDOUBLE(ISTART) points to the first doubling of DX that happens after JSTART
C  JDOUBLE(ISTOP) points to the last doubling of DX that happens before JSTOP
      
C    integration loop
      do i = istart,istop
         j1 = max(jstart,jdouble(i))+1
         j2 = min(jstop,jdouble(i+1))
         do j = j1,j2
            f3 = ucentr(j)+cntfug(j)-ecmn
            t3 = 2d0 * t2-t1+h1 * f2 * s2
            s3 = t3/(1d0-h2 * f3)
            reg(j) = s3      
      
            t1 = t2
            t2 = t3
            f0 = f1
            f1 = f2
            f2 = f3
            s0 = s1
            s1 = s2
            s2 = s3
         end do
         dx = 2d0 * dx
         h1 = dx * dx
         h2 = h1/12d0
         t2 = s3 * (1d0-h2 * f3)
         t1 = s0 * (1d0-h2 * f0)
      end do
      return
      end

